1. Field of the Invention
The present invention relates to computer systems for performing sparse matrix computations. More particularly, the present invention relates to a method and an apparatus that uses a hybrid approach to efficiently perform a cmod operation in solving a system of linear algebraic equations involving a sparse coefficient matrix.
2. Related Art
The solution of large sparse symmetric linear systems of equations constitutes the primary computational cost in numerous applications, such as finite-element design, linear programming, circuit simulation and semiconductor device modeling. Efficient solution of such systems has long been the subject of research, and considerable progress has been made in developing efficient algorithms to this end. A direct solution technique known as xe2x80x9cCholesky Factorizationxe2x80x9d is the most widely used approach to solve such a system. Under Cholesky factorization, the complete solution sequence requires many stages, including matrix reordering, symbolic factorization, numerical factorization and triangular solution. Of these stages, numerical factorization is typically the most computationally expensive.
One method of performing numerical factorization is based on a right-looking supernode-supernode method described in xe2x80x9cParallel Algorithms for Sparse Linear Systemsxe2x80x9d by Michael T. Heath, Esmond Ng and Barry W. Peyton, in xe2x80x9cParallel Algorithms for Matrix Computationsxe2x80x9d by Gallivan, et al. (Editors), SIAM (1994) (referred to as HNP). In a sparse matrix, a supernode is a set of contiguous columns that have essentially the same sparsity structure. Supernodes can be used to organize the numerical factorization stage around matrix-vector (supernode-column) and matrix-matrix (supernode-supernode) primitive operations leading to a substantial performance improvement arising from more efficient use of the processor caches and pipelining units.
The numerical factorization step involves two fundamental sub-tasks:
(1) cdiv(j): division of column j of factor by a scalar; and
(2) cmod(j,k): modification of column j by column k, k less than j.
These sub-tasks can be organized around supernodes. For example, cdiv(j) can be organized as an internal factorization/update of supernode j, and cmod(j,k) can be organized as a modification of supernode j by supernode k, k less than j.
Typically, the second sub-task is where the bulk of the computational cost is incurred. In order to increase computational efficiency, the cmod(j,k) operation can be divided into three steps:
(a) computation of the update and accumulation into a temporary array;
(b) carrying out the non-zero index matching between the first columns of the source and destination supernodes and computing relative indices; and
(c) scattering updates from the temporary vector into the target destination supernode.
By dividing the cmod(j,k) operation in this way, it is possible to apply techniques used in dense matrix operations in the step (a). Note that step (a) is where the dominant amount of time is spent. In the discussion that follows, we refer the step (a) as the xe2x80x9clocal dense cmod operationxe2x80x9d. The local dense cmod operation involves computing a trapezoidal-shaped dense update that can be represented as a combination of a dense rank-k update and a matrix multiplication.
Library routines can be used to speed up the cmod computation. These library routines are typically written in assembly language and are hand-tuned for a specific machine architecture. For example, on the current generation of UltraSparc-II-based machines, the Sun Performance Library (see http://www.sun.com/workshop/performance/wp-perflib/) provides SPARC assembly-language implementations of BLAS1, BLAS2 and BLAS3 routines. These hand-tuned assembly language implementations can yield performance close to the theoretical peak of the underlying processor.
For example, portions of the cmod operation can be efficiently performed by invoking the BLAS3 xe2x80x9cdgemmxe2x80x9d matrix multiplication code from Sun Performance library. Unfortunately, invoking the BLAS3 dgemm matrix multiplication code requires supernodes to be copied into and out of temporary storage because of incompatibility between data-structures used by a typical sparse solver and those expected by the dgemm API. This copying can add a significant overhead. Consequently, using the BLAS3 xe2x80x9cdgemmxe2x80x9d matrix multiplication code only makes sense for supernodes above a certain size. Otherwise, the performance gains from using the BLAS3 xe2x80x9cdgemmxe2x80x9d library code are cancelled out by the additional overhead involved in copying.
Hence, what is needed is a system that performs the cmod operation using library routines in cases where the performance gains from using the library routines exceed the computational overhead required to use the library routines.
Another difficulty in attaining high performance in numerical factorization is due to the fact that supernodes can have varying shapes, sizes, and sparsity patterns. These varying shapes, sizes and sparsity patterns can greatly influence computational performance. In order to optimize computational performance for the cmod operation, the supernodes of varying shapes and sizes must be divided into smaller sub-units for computation so as to balance computational operations with memory references in a way that is tuned for the particular machine architecture on which the computation is to be run.
One embodiment of the present invention provides a system for efficiently perform a cmod operation in solving a system of linear algebraic equations involving a sparse coefficient matrix. Note that CMOD and CDIV are used interchangeably with cmod and cdiv. Both cmod and cdiv are described below in conjunction with FIG. 3. The system operates by identifying supernodes in the sparse matrix, wherein each supernode comprises a set of contiguous columns having a substantially similar pattern of non-zero elements. In solving the equation, the system performs a CMOD operation between a source supernode and a destination supernode. As part of this CMOD operation, the system determines a subset of the source supernode that will be used in creating an update for the destination supernode. The system partitions the subset into a plurality of tiles, each tile being a rectangular shape of fixed dimensions chosen so as to substantially optimize a computational performance of the CMOD operation on a particular computer architecture. The system computes a corresponding portion of the update for the destination supernode using a CMOD function for each tile in the plurality of tiles. This computation may involve using computer code that is specifically tailored to each tile""s dimensions.
In one embodiment of the present invention, the system pre-computes a computational performance for a plurality of different tile sizes.
In one embodiment of the present invention, while partitioning the subset of the source supernode into the plurality of tiles, the system covers as much of the subset as possible using a first tile shape having a best computational performance. If there is a remainder after the first tile shape has been used, the system covers as much of the remainder as possible using a second tile shape with a second best computational performance.
Thus, based on the size of the supernode being processed, the system selects a substantially optimal tiling for computing the update.